---
title: A new automated chilled adult release system for the aerial distribution of
  sterile male tsetse flies
author: Caroline Mirieri1,4, Gratian N. Mutika1, Jimmy Bruno2  Momar Talla Seck32,
  Baba Sall43, Andrew G. Parker1, Monique M. van Oers54, Marc J.B. Vreysen1, Jeremy
  Bouyer1,65, Adly M. M. Abd-Alla1*
date: "29/06/2020"
output: word_document
---

```{r}
setwd("C:/BSI 2020/Submitted_final/raw_data")

library(ggplot2)
library(gcookbook)
library(ggfortify)
library(datasets)
library(MASS)
library(survival)
library(rmarkdown)
library(knitr)
library(coxme)
library(lme4)
library(nlme)

```


Figure 2a :Impact of motor speed (rpm) on the error rate in counting the releases males using the BSI machine

```{r}
tab2 <- read.csv("Fig 2a.csv")
head(tab2)
boxplot(tab2$error_rate~tab2$Speed,xlab = "Speed (rpm)", ylab = "Error_rate (%)")
fm2 <- lm (error_rate ~ Speed +(1|Replicate), data = tab2)
summary(fm2)
```

Figure 2b : Impact of flies density on the error rate in counting the releases males using the BSI machine at low speed (10 rpm)

```{r} 
tab2b <- read.csv("Fig 2b.csv")
head(tab2b)
boxplot(tab2b$error_rate~tab2b$Densities,xlab = "Initial no. of flies loaded into the machine", ylab = "Error_rate(%)")
fm2b <- lm (error_rate ~Densities+(1|Replicate), data = tab2b)
summary(fm2b)
```

Figure 3c:  Flowablity of BSI release machine: Prediction of manual counts from estimated counts 

```{r}
tab3b <- read.csv("Fig 3c.csv")
head(tab3b)


lm1 <- lme(nbcountedmanual ~ nbestim,
           random = ~ 1 | rep,
           data = tab3b,
           method = "ML")
summary(lm1)

plot((nbcountedmanual) ~ fitted(lm1), data = tab3b, pch = 1, cex = 2)
abline(lm((nbcountedmanual) ~ fitted(lm1), data = tab3b), col = "red")
cor.test(tab3b$nbcountedmanual,fitted(lm1))
summary(lm(fitted(lm1)~ (tab3b$nbcountedmanual)))$r.squared
```

###Impact of replicates on release consistency

```{r}
tab2 <- read.csv("flowability_replicates.csv")
head(tab2)
boxplot(tab2$no.per.min~tab2$Treatment,xlab = "T", ylab = " consistency of release rate")
fm2 <- lm (no.per.min ~ Treatment +(1|Rep), data = tab2)
summary(fm2)

```

Figure 4. Impact of release process treatment on male flight propensity 

```{r}
tab4 = read.csv("Fig 4.csv")
head(tab4)
tab4$rate <- tab4$out / (tab4$in. + tab4$out)
boxplot(tab4$rate ~ tab4$Treatment,xlab = "Treatment", ylab = "Flight propensity (%)")
tab4$Treatment <- relevel(tab4$Treatment, ref = "2-Irrd")
fm4 <- glmer(cbind(out,in.) ~ Treatment +(1|Replicate), family = binomial, data = tab4)
summary(fm4)

```

Figure 5. Impact of release process treatment on male Relative Mating Index(RMI)

```{r}

tab5 <- read.csv("Fig 5.csv")
head(tab5)
tab5$rate <- tab5$Treated / (tab5$Treated + tab5$Control)
boxplot(tab5$rate ~ tab5$Treatment,xlab = "Treatment", ylab = "Relative Mating Index (RMI)")
tab5$Treatment <- relevel(tab5$Treatment, ref = "2-Irrd")
fm5 <- glmer(cbind(Treated,Control) ~ Treatment +(1|Replicate), family = binomial, data = tab5)
summary(fm5)
```

Figure 6a. Impact of release process treatment on male premating period

```{r}
tab6a <- read.csv("Fig 6a.csv")
head(tab6a)
boxplot(tab6a$Premating_duration~tab6a$Treatment,xlab = "Treatment", ylab = "Premating duration (minutes)")
tab6a$Treatment <- relevel(tab6a$Treatment, ref = "2-Irrd")
fm6a <- lm (Premating_duration~Treatment +(1|Replicate), data = tab6a)
summary(fm6a)
```

Figure 6b. Impact of release process treatment on male mating duration. 

```{r}
tab6b <- read.csv("Fig 6b.csv")
head(tab6b)
boxplot(tab6b$Mating_duration~tab6b$Treatment,xlab = "Treatment", ylab = "Mating duration (minutes)")
tab6b$Treatment <- relevel(tab6b$Treatment, ref = "2-Irrd")
fm6b <- lm (Mating_duration~ Treatment +(1|Replicate), data = tab6b)
summary(fm6b)
```

Figure 7a: Impact of release process on the insemination rate of mated females

```{r}
tab7a <- read.csv("Fig 7a.csv")
head(tab7a)
boxplot(tab7a$MSV~tab7a$Treatment,xlab = "Treatment", ylab = "Mean Spermathecal Value")
tab7a$Treatment <- relevel(tab7a$Treatment, ref = "2-Irrd")
fm7a <- lmer(MSV ~ Treatment +(1|Replicate), data = tab7a)
summary(fm7a)
```

Figure7b:Distribution of spermathecal fill within treatments

```{r}
tab <- read.csv("spermfill.csv")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("spermfill-Control.csv")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("spermfill-Irrd.csv")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("spermfill-120m.csv")
summary(tab)
head(tab)
chisq.test(tab[,2:6])

tab <- read.csv("spermfill-5m.csv")
summary(tab)
head(tab)
chisq.test(tab[,2:6])
```

Figure 8A:Survival of released Glossina palpalis gambiensis under starvation conditions- KMA

```{r}
data_8A=read.csv("Fig 8a.csv")

#data_8A
str(data_8A)
head(data_8A)
data_8A=na.omit(data_8A)
names(data_8A)


#Kaplan Meier Analysis (km)
km <- with(data_8A, Surv(Day, Status))
head(km,15)
km_trt_fit <- survfit(Surv(Day, Status) ~ Treatment, data=data_8A)
plot(survfit(Surv(Day,Status)~Treatment,data=data_8A), xlab = "Day", ylab = "Survival rate", col=c('black','red','blue','green','orange'), lwd=2, xlim =c(0, 16)) 
legend('topright', text.font = , cex=0.8, c("1-Control","2-Irrd", "3-5m",  "4-60m","5-120m"), col=c('black','red','blue', 'green','orange'), lty = 1, lwd=2, box.lty = 1)

### survival: median Kaplan-Meier estimator 
survfit(Surv(Day,Status)~Treatment,data=data_8A)
data_8A$Treatment <- relevel(data_8A$Treatment, ref = "2-Irrd")
cox.surv <- coxph(Surv(Day, Status) ~ Treatment, data = data_8A, method ="exact")
summary(cox.surv)

```


Figure 8b:Survival of released Glossina palpalis gambiensis under feeding conditions-KMA

```{r}
data_8B=read.csv("Fig 8b.csv")

#data_8B
str(data_8B)
head(data_8B)
data_8B=na.omit(data_8B)
names(data_8B)

#Kaplan Meier Analysis (km)
km <- with(data_8B, Surv(Day, Status))
head(km,80)
km_trt_fit <- survfit(Surv(Day, Status) ~ Treatment, data=data_8B)
plot(survfit(Surv(Day,Status)~Treatment,data=data_8B), xlab = "Day", ylab = "Survival rate", col=c('black','red','blue','green','orange'), lwd=2, xlim =c(0, 65)) 
legend('topright', text.font = , cex=0.8, c("1-Control","2-Irrd", "3-5m",  "4-60m","5-120m"), col=c('black','red','blue', 'green','orange'), lty = 1, lwd=2, box.lty = 1)

### survival: median Kaplan-Meier estimator 
survfit(Surv(Day,Status)~Treatment,data=data_8B)
data_8B$Treatment <- relevel(data_8B$Treatment, ref = "2-Irrd")
cox.surv <- coxph(Surv(Day, Status) ~ Treatment, data = data_8B, method ="exact")
summary(cox.surv)
```





